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Abstract 

We report results of quantum Monte Carlo calculations of the ground state of dilute Fermi gases 
with attractive short range two-body interactions. The strength of the interaction is varied to 
study different pairing regimes which are characterized by the product of the s-wave scattering 
length and the Fermi wave vector, akp. We report results for the ground state energy, the pairing 
gap A and the quasiparticle spectrum. In the weak coupling regime, 1/akp < —1, we obtain BCS 
superfluid and the energy gap A is much smaller than the Fermi gas energy Epc- When a > 0, the 
interaction is strong enough to form bound molecules with energy Emol- For l/akp ^ 0.5 we find 
that weakly interacting composite bosons are formed in the superfluid gas with A and gas energy 
per particle approaching [£^^0^1/2. In this region we seem to have Bose-Einstein condensation 
(EEC) of molecules. The behavior of the energy and the gap in the BCS to EEC transition region, 
—0.5 < 1/akp < 0.5 is discussed. 

PACS numbers: 03.75.Ss, 05.30.Fk, 21.65. +f 
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I. INTRODUCTION 



How pairing evolves from the bare interaction has been a major question in condensed 
matter physics, and the study of pairing in relation to the phenomenon of superfluidity and 
superconductivity can be traced back to Cooper et al. 1]. Pairing lies at the core of several 
quantum many-body problems, and it is also believed to influence the evolution of neutron 
stars P]. Here we report results of quantum Monte Carlo calculations of a superfluid Fermi 
gas with short range two-body interactions. The strength of the interaction is varied to 
study different regimes of pairing. 

The evolution of pairing with the strength of the interaction has been discussed in the 
literature 0], 0]. In the regime where the interaction is weak and attractive, a gas of 
fermions has a superconducting instability at low temperatures, and a gas of Cooper pairs is 
formed. The typical coherence length is larger than the interparticle spacing tq {AirrQp = 3 
with p the number density) and the bound pairs overlap. In contrast, in the strong-coupling 
limit the coherence length is small, and the bound pairs can be treated as well seperated 
Bose molecules. One then expects the molecules to undergo Bose-Einstein condensation 
(BEC) into a single quantum state with zero momentum. 

I — > 1—1 

The Bardeen-Cooper-Schrieffer(BCS) theory 0] and Gorkov equations j;5] have been used 

to estimate gaps in superfluid gases. However, their predictions differ by more than a factor 
of two and they may be qualitatively valid only in the weakly interacting regime. Here we 
use first principle quantum Monte Carlo methods to study the entire region ranging from 
free fermions to the tightly bound Bose molecules. 

Dilute Fermi gases of ^''K, ^Li, ^H for example, can now be studied in the laboratory using 
magnetic and optical trapping and ingenious cooling methods Q], These are dilute Fermi 
systems, in contrast to dense atomic liquid ^He or a solution of ^He in superfluid ^He. Within 
the last few years temperatures T <^Tf have been achieved, where Tp = is the Fermi 
kinetic energy and kp is the Fermi wave vector. At such a low temperature, the fermionic 
nature of the quantum statistics becomes evident in the measurement of the density profile of 
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the trapped gas. At even lower temperatures the transition to the superfluid Cooper-paired 
state is expected. However, the temperature Tc of this transition can be much lower than Tp 
and conclusive evidence of superfluidity is still to be seen. In order to have the transition at 
an achievable temperature, the experimentalists rely on the Feshbach resonance technique 
to produce strong interaction between the fermionic atoms. 



When the range of the interatomic interaction is smaller than all the length scales in 
the system, the details of the interaction are believed to be unnecessary and the scattering 
length a is sufficient to characterize it. Near the resonance, the magnitude of the scattering 
length a becomes much larger than Tq and the system enters the strong-coupling regime. 
The value akp ~ —7.4 has been achieved by O'Hara et al. f^] and the limit {akp —oo) is 
now approached in the laboratory 
atoms was reported by Regal et al 



(Ref. (ll|) 



. Recently, creation of bosonic molecules from ^''K 
10|, and pairing in the 1/akp ~ regime was observed 



A few words are in order regarding the language of s-wave scattering. For a noninteracting 
system at zero temperature, the only length scale is 1/kp. We can use the dimensionless 
quantity akp to describe a dilute gas having interparticle spacing tq much greater than 
the interaction range. We often use 1/akp because akp changes discontinuously from —oo 
to +00 when a bound state is formed at 1/akp = 0. For attractive interactions 1/akp can 
change from large negative values (weakly interacting limit) to large positive values (strongly 
interacting limit). As discussed in section m the radius of the bound molecule provides 
another length scale in the strongly interacting regime. Some physical examples of the limits 
of 1/akp are: 1) electrons in superconductors have 1/akp large and negative; 2) neutron 
matter has 1/akp small and negative; and 3) cold deuterium atoms have large positive 1/akp. 
In the last case, molecular bound states smaller than the average interparticle distance tq are 
possible. On the other hand, superfluid ^He is not describable in terms of akp, because the 
interaction range is greater than ro, and the paired state does not have s-wave symmetry. 

In the limit of zero energy for the colliding pair, the two-body scattering cross section 
a is given by 47ra^. When \a\ -C ro, the interatomic collisions in the gas are similar to 
those in vacuum, and the mean free path is approximately given hj i = However, this 
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approximation is meaningful only when \a\ <^ tq and i > tq. When |a| is ^ ro the two-body 
collisions in the gas are strongly influenced by the presence of other particles, and their cross 
section in the gas is much smaller than in vacuum. 

For a Fermi gas at low density, an expansion of the energy in terms of akp is possible. 
For spin 1/2 Fermi gases it is known to be 



10 4 

1 + —{akp) + —(11 - 21n2)(aM' + Oiakp)^ + ... 
9tc 21 



(1) 



where Epc = f ~2nf" ~ f-^^ ground state energy per particle of the noninteracting 

Fermi gas. In the akp —>■ — oo limit, theoretical estimates of 0.326 and 0.568 Epc were 
reported [l^, More recently, the authors |l^ predicted Eq = (0.44 ± 0.01)i?irG using 
quantum Monte Carlo methods. In this paper we continue that study of the properties of 
cold dilute spin 1/2 fermion gas and extend it to all the regimes of 1/ akp a.s a. first step for 
understanding the superfluidity and the bosonization of dilute Fermi gases. 



The model considered in this study consists of A fermions contained in a box with periodic 
conditions on its boundaries. It is not polarized so that half of the spins point up and the 
other half down. Typically A is varied from 10 to 20 to estimate properties of uniform gas 
in the A ^ oo thermodynamic limit. In some cases larger values of A are used. Fermions 
of the same spin do not feel the effects of interaction because it is of short range and Pauli 
exclusion predominates. The fermions of different spins interact via a central potential v{r) 
with the following properties: 1) It is attractive with very short range as we assume the 
dilute limit, 2) The details of the potential do not matter, in principle we can think of it as 
an attractive delta function potential, and 3) The potential can be adjusted such that we 
can sweep through different regimes of akf. 

From the considerations mentioned above, a cosh potential of the form 

v{r) = -vo -J- — - , (2) 

m cosh [fir) 

can be used. The strength of potential (vq) is adjusted to obtain the desired value of akp. 
We can also take appropriate values of /i such that the effective range of the potential -Re// 
is much smaller than the interparticle distance Tq. When f o = 1 this potential has a = ±oo 



and Reff = 2//i. In most calculations we have used /iro = 12. For the a — oo case we 
also tested the ^r^ — > oo limit using /xro = 24 |l6| . 

Results of simple lowest order constraint variational (LOCV) calculations are reported in 



sectio n ITT| The LOCV method was first used to study neutron matter 
et al. 
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. Recently, Cowell 



18| have used it to study cold Bose gases in the unstable a > regime. It provides 
a surprisingly good estimate of the ground state energy. Here we use it to study the effect 
of the difference between the cosh (yuro = 12) and delta-function potentials on the energy 
of dilute gases. The difference becomes significant when 1/akp — > cxd, and the radius of the 
molecule approaches 1/ ^. LOCV is also used to estimate the energy of the unstable state of 
the Fermi gas for a > 0. The stability of dilute gases is discussed in the LOCV section ITTl 

One of the limitations of LOCV is that it can not be used to calculate the pairing energy 
gap A or the other superfluid properties of Fermi gases. The quantum Monte Carlo methods 
used in Ref. |1h] and this work to study superfluid gases are described in section UTTl and the 
results for the energy, pairing gap and the quasiparticle spectrum are presented in section 
II VI over the range akp = —1 to =Foo to +0.5. Conclusions are given in the last section IVl 



II. LOWEST ORDER CONSTRAINT VARIATIONAL CALCULATIONS 



In the lowest order constraint variational (LOCV) method the ground state of the Hamil- 
tonian 

^ = -|^Ev; + E-K'), (3) 

where the unprimed index i denotes spin up particle, primed index j' denotes spin down 
particle, and p can be any particle, is approximated by the Jastrow-Slater wave function 

i^w=n/K')i*5) ' (4) 

where \^s) is the ground state of noninteracting fermions. In the present case \^s) is a 
product of two Slater determinants, the first corresponding to the spin up fermions and the 
second corresponding to the spin down fermions. The interaction effects are represented by 
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the Jastrow function H fi'f'ij'), where f{rij>) denotes the pair correlation function. We often 
use fij> to denote f{rij'). fiji = 1 means no correlation between the pair ij' and fij> ^ 1 for 
correlated pairs. In variational calculations, the function /(r) is determined by minimizing 
the expectation value of the Hamiltonian 



-t^i'^si n fap' VI n M'^s) + E($5i n fa^' n M'^s) 
^ '~ {^s\ n fap' n M^s) ■ ^ ' 

The assumption behind LOCV is that the energy is most sensitive to the correlations 
of short (less than tq) range. We impose a constraint on the range of /(r) to assure that 
the correlations are mostly among the closest pairs, and keep only the pair terms in the 
cluster expansion of the energy expectation value. The healing distance d is the range of 
/(r) defined such that /(r > d) = 1 and ^^^\r=d = 0. In LOCV, d is chosen such that 
on average there is only one other particle within the distance d of any particle. Effects of 
deviations from this average are assumed to cancel. 

Euler-Lagrange minimization of the energy expectation value gives a Schrodinger-like 
equation for f{r<d) 

--V'f{r) + v{r)f{r) = \f{r) , (6) 
m 

The constraint used to determine the healing distance is 

P-j'nr)dh = l, (7) 

and the A is chosen such that = 0. In the equations (jHI) and ((Tj) we do not have 

exchange contributions because the range of the interaction is short and fermions of same 
spin do not interact. When the equations © and dTj) are simultaneously solved, the energy 
per particle is given by 

Elocv = EpG + — . (8) 

The results obtained for the ground state energy of spin 1/2 Fermi gas with the cosh and 
delta-function potentials are shown in Fig. ^ 
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When l/akp < the ro is the only length scale in the gas, and the results obtained 
with the cosh potential with firo = 12 are indistinguishable from those given by the delta- 
function potential. In contrast, when 1/akp > 0, we have a molecular bound state whose 
radius provides another length scale. At large positive values of 1/akp there are differences 
between results of the present cosh and delta-function potentials due to the rms radius, 
Rrms of the molecule becoming comparable to the range of the present cosh potential. For 
example, at l/akp = 2 we get fiRrms = 2.3 with the present choice of /i. In principle, 
we can continue to approximate the delta-function interaction with the cosh potential by 
further increasing /i, and working in the ixRrms ~^ oo limit. However, all of the present 
computations are with (ivq = 12. 

Fig^also shows the presumably exact results obtained with the cosh potential with the 
GFMC method described in the next section. The LOCV energies appear to be surprisingly 
accurate. However, it should be realized that a part of the accuracy of LOCV is due to 
a cancellation of errors, and not due to the quality of the Jastrow- Slater variational wave 
function (Eq. HJ. In fact, the variational energy upper bound obtained with that wave 
function for l/akp = is = (0.62 ± 0.01)^^, significantly above GFMC result of (0.44 ± 
0.01)Efg- The LOCV energy of OAGEfg is below the Jastrow-Slater variational upper 
bound because it is calculated approximately keeping only two-body cluster contributions. 
However, when the contributions of > 3-body clusters become important we can expect that 
the approximations in the Jastrow-Slater wave function would also become important, and 
the true energy will be below the Jastrow-Slater upper bound. 

The ground state energies obtained with the conventional BCS (variational) method are 
also shown in Fig. ^ In the weakly interacting limit, 1/akp —>■ — oo, the BCS energy is too 
large since it does not have the correct low density limit given by Eq. ^ On the other hand, 
in the strongly interacting limit, 1/akp —>■ +oo, the BCS energy is very close to the exact 
result (GFMC) presumably because in this limit we have complete pairing of the fermions 
into Bose molecules. LOCV is less accurate than the conventional BCS method in strong 
coupling region. 
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The LOCV pair correlation functions are shown in Fig. |21 The heahng distance d ^ ro in 
the weakly interacting region {1/akp << 0), and as we increase the strength of the potential, 
/(r) becomes more and more peaked at the origin, and d becomes smaller than rg. In fact 
for 1/akp » 0, the boundary condition at d has less impact on Epocv and A of the Eq. El 
becomes close to the molecular binding energy Emoi such that Epocv = EpG + + ^E. 

is the term that predominates in this limit. 6E is small {\SE\ < Epc) and negative so 
that Epocv > Emoif^- 

When a > we can obtain another solution of the LOCV equation with a node at r < d. 

1 — I 

This solution was discussed by Cowell et al. for cold Bose gases, and at small values 
of akp it gives results in agreement with the low density expansion, (Eq^. The first term 
iiak,) is conectly reproduced by LOCV, but the higher order terms are apppximate. 
In the limit a — > oo we have the condition kdta.n{kd) = —1 discussed in Ref. The 
solution with one node is kd = 2.7983 and it gives E/N = Epc + f = Epc + ~ 
3.92EpG- Results obtained with the delta-function potential, including this unstable region 
are shown in Fig. El Those corresponding to the nodeless solution of the LOCV equation 
are represented by full line, while the dashed line corresponds to the solution with a node. 

The state of the gas having a node in the pair correlation function /(r) is unstable because 
it has energy > Epc, while that with nodeless /(r) has lower energy < Epc (see Fig. Oj). 
However, it can have a relatively long life time because energy conservation prevents two 
atoms to make the transition to the lower energy state. At least three atoms are needed, 
which hinders the transition at low densities. Most of the observed BEC of Bose atoms are 
in such unstable states in which the /(r) has nodes at small r. 

The E{akp) shown by the solid line in FigjHl corresponds to the stable ground state 
of the model Hamiltonian with the delta-function interaction. In principle, this state can 
be exactly calculated by the quantum Monte Carlo method described in the next section. 
However, when the range of the interaction is finite, as for the cosh model, the system can 
collapse to a tightly bound state at large density. This instability can be easily seen in the 
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Hartree mean field approximation in which 



Emf{p) = Efg{p) + ^Iv 



h = J v{r)d'r , (9) 

where 1^ (< 0) is the volume integral of the interaction. At large enough p the interaction 
energy becomes larger than Epc leading to a tightly bound state. 

Consider for example a simple square well potential of range R such that v{r < R) = — Vq 

fe 2 2 

and v{r > R) = 0. Let this potential correspond to a = oo. This means Vq = ^^^2 and, 
/„ = -^VoR^ = -^R. Then 



3 n'kl L Stt 

^mf{Kf 



5 2m 



1 - ^fl*. 



(10) 



The collapse occurs at values of kp > |^-^, and can be pushed to higher densities by reducing 
R, or equivalently increasing /i in the case of the cosh potential. In the present studies, we 
ignore this collapsed state; assuming that it occurs at too large a density to influence the 
dilute gas properties. 



III. GREEN'S FUNCTION MONTE CARLO CALCULATIONS 

Green's function Monte Carlo (GFMC) [2^ is a powerful method for calculating the 
ground state properties of many-body quantum systems. It can be used to calculate the 
ground state properties of Bose systems with controllable statistical errors without approx- 
imation. For the fermion systems, however, we have to deal with the sign problem posed by 
the anti-symmetry of the wave function as discussed below. We begin with a brief overview 
of the GFMC method. 

Let be the eigenstates of 7i with eigenvalues Ei. The trial variational wave function 
which provides an approximation to the ground state \l/o, can be expanded as 

^y = E«^^- (11) 
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In GFMC we project out \I/o from by evolution in imaginary time 

^(r ^ oo) = lim e-(^-^^)"^i/ = lim V a^e-^^'-^^^'^i 

T— >00 T^OO 

i 

— ^ aoe-^^^-^^^'^o , (12) 

where we have shifted the origin of energy to Et ~ Eq to control the norm of \l/(r ^ oo). 
In practice, the time evolution is carried out in n small steps 



e 



-{H-Et) 



r = ^^-{n-Er)i^r ^ Ar = r/n, (13) 



and Et is tuned to keep (\l/(r)|^(r)) constant. The tuned Et provides the growth estimate 
of the true Eq. An alternative method for calculating the ground state energy, often with 
smaller statistical error, is given by the mixed estimate (see Fig. 

= {^>ymr ^oo)) = i^yMr oc)) = ■ 

In general, the time evolution operator or propagator is not known for arbitrary large 
value of r except for few simple systems. However, we can obtain small time prop- 
agator with controllable errors for any Hamiltonian with static potentials that depend 
only on the positions of the particles denoted by a 3A^-dimensional configuration vector 
R = {ri, r2, . . . ; r'^^, r'2, . . .}. This is the motivation to write the time evolution as a product 
of many short time operators (Eq. fT!?|l . We define the Green's function 

G'(R,R') = (R|e-(^-^^)^"|R') . (15) 

The propagation equation becomes 

^(R, r + Ar) = J dR'G(R, R')*(R', r) . (16) 

The primitive approximation to this Green's function is 

G(R, R') ^ e-(^W-^-)^Go(R, R')e-(^(^')-^-)^ , (17) 

where ^(R) = I] '^ij'ij') and Go(R5 R') is the Green's function for A free particles 



Go(R, R') 



m 



^ 3 



2-Kh'AT 



A 



2 

e 



(18) 
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This approximation has errors of order Ar^. The total error after n time steps is of the 
order ~ nAr'^ = The corrections to this expression can be sampled to make an exact 
algorithm. Here we use the more common method and make this error as small as we 
want, by increasing the number of steps n. In practice, this error is made smaller than the 
statistical sampling errors of the Monte Carlo integration. 

A naive quantum Monte Carlo algorithm could start with A^^ configuration vectors Rj 
sampled from |\l/y|. These provide the approximate representation 



where wi = 1 or —1 depending on the sign of The accuracy of this representation 
increases with the number of samples A^^. Inserting Eq. into Eq. and using the 
short time approximation, gives \E'(R, Ar) as a sum of normalized gaussians times weight 
factors containing the product of the original Wi and the exponentials in short time Green's 
function (Eq. fTTj) . Sampling a position from each of the Ns gaussians gives a representation 
of \1/(R, Ar) as a sum of delta functions times weight factors with signs. This process is 
repeated n times to obtain \1/(R, r = nAr). During the evolution, large magnitude weight 
factors are converted into multiple copies while small factors are sampled and kept with 
unit magnitude new weight with a probability proportional to the magnitude of the old 
weight. The random walk of the weighted 5-function samples representing the propagation 
of \1/(R, t), therefore consists of diffusing and branching and the number of samples at each 
time step can vary. 

This algorithm suffers from the fermion sign problem. For A^s samples Ri, 1 < « < A^s, 
and weights Wi, the denominator of a matrix element such as the mixed energy will be the 
sum J2i^Ji w^i^y(Rj)- Each Wi carries the sign of the initial sample from and if the path 
of the sample i has crossed nodes of \I'v' odd number of times, the contribution to the sum 
will be negative. For large times the contribution of these negative paths almost completely 
cancel the contribution of positive paths that have not crossed nodes or crossed an even 
number of times. The signal dies out exponentially compared to the statistical noise. The 
numerator suffers from the same problem. 




(19) 



i=l 



11 



The fixed node[21[ approximation deals with the fermion sign problem by restricting the 
path so that crossings of the nodal surface are not allowed. When this constraint is imposed 
with the nodal surface of the exact fermion ground state, \E'(R, r) converges to that state. 
Imposing the nodal surface from an antisymmetric trial function gives an upper bound 

lim {n)mix,r > Eo . (20) 

We impose the fixed-node constraint with the nodes of our trial function \E'y(R). 



Importance sampling is used to control the fluctuations of the weights. The propagation 
equation is modified by multiplying by a positive importance function. Since we are using 
the fixed node approximation, the paths have zero probability of crossing the nodes, and we 
can take the importance function to be the absolute value of The propagation equation 
now becomes 

[|vI/^(R)|M/(R,r + Ar)] = J rfR'^^G(R, R') [|v&y(R')l r)] . (21) 



The short time approximation for the importance sampled Green's function is 

^GofR,R' + ^vin|MR^-^J IMR)I^^(R,R' 



2m ; [|^v'(R')|G'o(R,R' + ^Vln|^v'(R) 

^ Go (^R,R' + ^Vln|vl>^(R)pj |e-(|[i^.(R)+i^.(R')]-i^T)A.| (22) 

where the local energy is 

^'^^^ - ^7(sr ^^^^ 

Since Go is still a normalized gaussian, the only changes to the naive algorithm are the 
sampling of the drifted Gaussian, and the new weight given by the terms in the braces. 
Notice that if is close to the ground state of Ti, El(R) will have less fluctuations than 
l^(R), and the branching of the walk is much reduced. Any paths that cross a node due to 
the short time approximation are eliminated. 

For Ng samples Rj, all with weight Wi = 1, at time r, the mixed energy becomes the 
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average of the local energy 

(H)™,. = . (24) 



' s 



Since the fixed node calculations give an upper bound to the ground state energy, our 
strategy (see Ref. is to choose a trial wave function with variable nodal surfaces and 

minimize the fixed node GFMC {'H)mix- 

The trial wave function \l/y(R) is now used in three different contexts; 1) as the initial 
guess of the ground state, 2) as the importance function in Eq. EH and 3) as the node 
restriction function. The nodes of the Jastrow-Slater wave function (Eq. equal those of 
noninteracting Fermi gas and can not be varied. So that wave function is not useful for 
present studies. 

From physical considerations, a better trial wave function must reflect the fact that the 
fermions with attractive interaction can form bound Cooper pairs in the ground state. And 
from mathematical considerations, the trial wave function must have variable nodal surface, 
which can be varied to minimize the flxed node GFMC energy. The BCS wave function is 
such a wave function. Commonly, we write 



\BCS)e = n(^k, + e*^k,aL,t«-k4)|0) , (25) 
p 



where |0) denotes the vacuum and u\^^ and are real positive numbers. However, this 
wave function does not correspond to a deflnite number of particles. In fact, expanding the 
wave function we can write 

\BCS)e = \Q) +e''P^|0) + e'^\P^f\Q) + e'^\P^f\Q) + ... . (26) 

where = E ^^kpt^-kpi P^^^ creation operator. The component that corresponds 

to A particles or M = A/2 pairs can be obtained by transforming 

\BCS)a = ^1^''''' \BCS)ede, 



= (pt)M|o) . (27) 
13 



This component can be written as an antisymmetrized product of the pair wave functions 



^BC5(R) = ^[0(riiO0(r22')-0(^MM')] , (28) 

p "^kp p 

where the number of up spin particles (M) is equal to the number of down spin particles 
(M'). The variational parameters ap are real positive numbers. The free fermion gas, Slater 
wave function is just a particular case of this wave function when ap ^ for |kp| < kp and 
= for |k„| > kp- 



We also consider systems having unpaired particles. In particular, we can have M pairs 
and 1 unpaired up or down spin particle. This generalization is necessary as the gap energy A 



is calculated from the odd-even staggering of the ground state energy 



16|. With 1 unpaired 



(t or I spin) particle in the state '?/'k^(r), with momentum k^^, the trial wave function is 



given by 
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^BCs(R) =^{[0(rnO...0(rMAf')]^k„(r)} . 



(29) 



The ground state is expected to have |ku| ^ |k^| in the weakly interacting regime and 
k„ — > in the strongly interacting regime. This wave function can be calculated as a 



determinant 



,0 



which makes the numerical calculations relatively simple. 



Quantum Monte Carlo calculations use a finite number of particles in a cubic periodic 
box of volume to simulate the infinite uniform system. The momentum vectors in this 
box are discrete 

271 

kp = —{rip^x + npyjj + Hpzz) , (30) 

and the system has a shell structure with closures occurring when the total number of 
particles = 2, 14, 38, 54, ... for spin-1/2 fermions. The shell number I is defined such that 
I = nl + nl + nl and = ^^/. 

In the present calculations, the pair wave function (/)(r) has the assumed form 

0(r) = ^(r) + ^1^"^''^ > (31) 
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/3(r) = p{r) + P{L - r) - 2/3(L/2) for r < L/2 







for r > L/2 



/3(r) = [1 + -fbr] [1 - e'"'''] 




Here = 4 is a cut off shell number. We assume that the contributions of shells with I > Ic 



range L/2. We further reduce the statistical fluctuations by using the Jastrow factor along 
with "^Bcs in the variational wave function: 



The Jastrow factor does not change the nodal structure. Thus, the average value of the 
estimated energy is independent of /(r), but the statistical error is reduced by using the 
/(r) from LOCV calculations. 

It is convenient to require that ^ = at r = 0. This is because the local energy has 
terms like which can have large fluctuations at the origin when 7^ at r = 0. The 
factor [1 — e"'^'"'] cuts off ^ dependence of /5 at 6r < i. The energies are not too sensitive 
to the parameter c, and its value is fixed at 10. In addition, 7 is chosen such that ^ = at 
r = 0; its value is 6 in the limit L 00. 

The variational parameters are {ao, ai, . . . , aj^} and b. We wish to find a set of these pa- 
rameters that minimize the fixed node GFMC estimate of energy. However, considering that 
we have to allow simultaneous variation of all the parameters, methods based on unguided 
variation become difficult, if not infeasible. Again, we rely on the GFMC procedure itself to 
optimize these parameters. Initial configurations are obtained with a random distribution of 
the parameters centered around a reasonable guess. Each of them is propagated according 
to the nodal constraints provided by their parameters with a single Et- The paths with the 
smallest {7i)mix acquire large amplitudes or weights as r — 00. The average among these 
paths gives an optimization over the initial random distribution. This process is repeated 
several times until convergence is achieved. 



to the pair wave function can be approximated by a spherically symmetric function j3{r) of 



^y(R)=n/(n.')*BC5(R) . 



(32) 



When we have an odd number of particles, the ground state momentum k„ (Eq. 0^ is 
an additional variational parameter. We minimize the fixed node GFMC energy of systems 
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with odd A by varying k^j. As discussed in the results section, the magnitude of k„ changes 
from fci? to as the interaction strength increases and we go from the weakly interacting 
BCS to strongly interacting BEC regime. The gap energy is obtained from the odd-even 
staggering of the total energy 

A(2M + 1) = E{2M + 1) - ^{E{2M) + E{2M + 2)) (33) 

In doing so, the effects of interaction among quasiparticles are neglected. 

Results for akp = oo are shown in Fig. |3| The energy per particle E/A and the gap A do 
not have a significant A dependence in this case. These results were reported in Ref. 
and results for other values of akp are presented in the next section. 



IV. RESULTS 



The values of the parameters, ao-4 and 6, of the BCS wave function are to be determined 
by minimizing the fixed node GFMC mixed energy for each value of akp and A. The 
minimum energy obtained is our estimate for the ground state energy of the system. The 
values of the parameters that minimize this energy are not very sensitive to A, the number of 
particles in the box. We find it sufficient to determine the optimum parameters at A= 10, 14 
and 20, and interpolate their values in the A = 10, 14 and A = 14, 20 ranges. The values of 
the parameters at these values of A are listed in Table HI 

At akp = —1 the lowest energies are obtained without any short range f3{r) and the 
optimum pair function has contributions only from the states with / < 3. This is consistant 
with the weak coupling BCS theory in which goes to zero when k — kp becomes large. 

When 1/akp > —1, lower energies are obtained with /5(r) 7^ 0. In most cases, the values 
of the parameters do not change significantly between A = 14 and 20. The values listed in 
Table m for A = U are used for U< A< 20. 

In the < 1/akp < 2 range, the optimum values of the parameters of "^bcs do not seem 
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to change significantly with the akp- We have not obtained any significant improvements to 
the energy from varying the parameters in the region l/akp > 0. In this region we retain the 
values found for 1/akp = and A = 14. Recall that only the nodal surfaces of the ground 
state wave function are constrained by those of "^bcs- The complete \l/y has an additional 
product of Jastrow pair correlation functions f{riji) which depends on akp, and the true 
ground state wave function changes continuously with akp- 

The magnitude of the momentum \^u of the unpaired particle in the ground state is also 
determined by minimizing the GFMC fixed node energy. The minimum values are listed in 
Table ITTl In the weak coupling limit, the BCS ground state for odd A has |ku| = Ikp^l. In 
the periodic box, the value of fc^ is 1 for 2 < A < 14, and 2 for 14 < A < 38 in units of 
{2'n/Lf. In the akp = to -3 range, the minimum values of k,, indicated by the 

weak coupling BCS theory. However, in the —10 to 3 range the fc^ is 1 for the entire range 
(11 to 19) of odd values of A considered. At akp = 2 the states with k"^ = and 1 are 
almost degenerate, and for l/akp > 0.5 the ground states of odd A systems have k„ = 0, 
as expected when the system consists of bound molecules condensed in the zero momentum 
state, and the unpaired particle also in the k^ = state. 

The calculated values of the ground state energy are shown in figures IHl and [7| The 
systems with l/akp < 1/3 seem to have E > 0, while those with l/akp > 1/3 can have 
E < 0. When 1/a > the two-body interaction is strong enough to bind two particles and 
form molecules with energy Emoi- The energy per particle, E/A of the superfiuid Fermi gas 
is compared with Emoi/2 in figure |H1 Within the computational errors E/A> Emoif^ (see 
Table UTT]) . however at l/akp > 0.5 we find that E/A is very close to Emoif^- This behavior 
also indicates that at these values of akp the system approaches that composed of Rose 
molecules forming a BEC. It has been argued that the interaction between these molecules 
is weakly repulsive, with a molecule-molecule scattering length given by a^m = 0.6a js^. 
In this case the E/A will always be greater than Emoi/'^i and the gas will have positive 
pressure, E/A increasing with the gas density oi kp- 

The pairing gaps calculated from the odd-even energy difference (Eq. are shown in 
figures 1^ and ITUl These gaps are not very sensitive to A, and they are compared with the 
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predictions of BCS jsl] and Gorkov estimates given by 



Abcs = -Tf e"/(''^'=^) 



e 

/Os 7/3 1 /9\ 1/3 

Aco..o. = (^) Tf e-/^'^'-^ = ^ (^) Ascs • (34) 

where the chemical potential is approximated by Tp as when l/akp << 0. At l/akp << 
the calculated gaps are in between these estimates, while at positive values of l/akp they 
approach —Emoi/'^ as expected for a gas of Bose molecules (see Fig. Uniand Table ITTT|l . 

Figures |H1 and ^2 also show the Emoi for a delta-function interaction in addition to those 
for the present cosh potential with /xro = 12. The two potentials give essentially the same 
results for 1/akp < 0.5, but at larger values the cosh potential is more attractive. The 
values of the rms radius Rrms of the molecule are listed in Table IIIII At large values of 
1/akp the fiRrms is not very large for the present choice of yurg = 12, and much larger values 
of /i should be used to approximate the delta-function interaction. 

The pressure (P = P^^^f^) and the adiabatic index (F = -p^) of the superfluid gas in 
the range —20 < akp < are shown in Fig. El For noninteracting Fermi gas {akp = 0), 
E/A = EpG and we have P = ^pEpc and F = |. In the limit akp — > we can use the low 
density expansion (Eq. [T)) to obtain 

2 5 

P{akp 0) ^ -pEpcil + —akp), 

3 Svr 

5 5 

T{akp -^0) ~ - + —akp. (35) 

3 9tt 

In the akp —>■ —oo limit, we have E/A = ^Epc, therefore 

2 

P{akp -oo) = -ipEpG 
5 

V{akp — oo) = - (36) 
o 

where ^ = 0.44 ± 0.01 according to the present calculations. 



The calculated value of P{akp) /{pEpo) is | at akp = and decreases monotonously 
to 1^ as akp — oo. However, The adiabatic index T{akp), is | for both akp = and 
akp = — oo, and has a minimum value of ~ 1.6 at akp ~ —1.3. 



18 



V. CONCLUSIONS 



The present work shows that accurate calculations of the pairing gaps and energies of 
superfluid Fermi gases are possible with the fixed node GFMC method. The unknown nodal 
surfaces can be determined variationally by minimizing the fixed node GFMC energy. This 
method gives the exact result in the l/akp —oo (Fermi gas) and l/akp +00 (BEG of 
molecules) limits for short range attractive interaction, and seems to overcome the fermion 
sign problem. An alternative method based on path integral Monte Garlo simulations is also 
being developed j^]. 

Our results are in qualitative agreement with the known BGS-BEG crossover model 
(see Leggett 0) where gap and chemical potential(/io) are calculated self consistently. 
The gap is determined as the minimum of the Bogoliubov quasiparticle energy E]^ = 



ek — /Uo)^ + l^kP) where ek = is the single particle excitation energy and is 
the gap parameter. Two limiting cases were considered in this referenced article. For 
l/akp —00, fiQ ^ Tp > 0, the minimum of E'k occurs at k = kp and the minimum quasi- 
particle energy A = A^^ = ^ Tp e^/^'^"-^''\ However, for 1/akp +00, /io ~ -E'mo//2 < 0, 
the minimum of is at /c = 0, and its value A = |/io| ~ \Emoi\/'^ because A'j^ ~ 0. The 
BGS-BEG crossover takes place when /xq = and this corresponds to akp positive and of 
the order 1. The odd-even staggering A(2M -|- 1) given by Eq. IHHl presumably equals the 
minimum quasiparticle energy in the limit M — 00. 

According to Leggett's description, in the weak BGS superfluids the ground state of 
systems with odd number of particles is expected to have momentum kp, while in the 
molecular liquid with BEG it is expected to have zero momentum. With this criteria the 
calculated values of /c^ (Table |n} suggest that the BGS to BEG transition occurs in the 
range —0.5 < 1/akp < 0.5. It appears to be a smooth transition or crossover. 

A recent experiment by Bartenstein et al. also seems to corroborate some of our find- 
ings. In fact, in their paper 2^ BGS-BEG crossover regime for ^Li is reported to be 
—0.5 ^ 1/akp ^ 0.5. In addition, in the unitary limit {akp = ±00) they measured 
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E/A = Q.Z2tfQEFG which includes within its range our result E/A = (0.44 ± Qm)EFG- 

We can notice that in the BCS regime A is much smaller than E/A^ while in the BEC 
regime A ~ I-E'IM ~ \Emoi\/'2- However, in the transition region A is significantly larger 
than \E\/A. 
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FIG. 1: Ground state energy per particle of dilute Fermi g ctSGS clS 8b function of akp- The full and 
dashed curves give the LOCV results for cosh (fiVQ = 12) and delta-function potentials, and the 
circles show the essentially exact results for the cosh potential obtained with the GFMC method 
described in section [TTTl The dotted and dash-dotted curves correspond to the conventional BCS 
results with cosh and delta-function potentials respectively. 
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FIG. 3: The LOCV E/A in units of EpQ vs akp for attractive delta-function potential. The dashed 
line corresponds to /(r) having one node, and the solid line shows the results with nodeless f{r). 
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FIG. 4: Mixed and growth estimates of the GFMC energy. The r — > oo asymptotic value is reached 
after ~ 5000 times steps. Each time step At is 1.2 10"^-^^. 
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FIG. 7: -E(^) for l/akp > 1/3. The results for akp = 0.5 have been multiplied by 0.5 for graphing 



27 




0.4 



0.8 1.2 

l/(akj 



1.6 2 



FIG. 8: E/A and Emoif^ for positive values of l/akp for the cosh (/xro = 12) potential and Emoif^ 
for the delta-function potential. 
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FIG. 11: Calculated values of ^qpmc{0'^f) (cosh, /xro = 12 potential) are compared with various 
estimates of ^{akp) and —Emoi/'^- The BCS and Gorkov estimates do not depend on the shape of 
the potential, while — £'moi/2 is shown for both cosh (solid line) and delta-function (dash double 
dot) potentials, ^bcs ^-nd ^Gorkov assume the chemical potential ~ Tp throughout the whole 
range of akF{see equations 151]). 
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FIG. 12: Pressure(P, dashed curve) and adiabatic index(r, continuous curve) in the BCS regime 
{akp < 0). In the dilute limit {akp — > 0) we have P/{pEpG) — > | and F — s- |. In the dense limit 
{akp — oo) we have P/{pEfg) ~^ 0.44| and F — > |. F has a minimum at akp ~ —1.3. 
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TABLE I: Optimum values of the variational parameters. 



akp A 


ao oti a2 OC3 a4 b 


-1 10 

14 
20 


1.00 0.05 NA 
1.00 1.00 0.010 NA 
1.00 1.00 0.104 0.024 NA 


-3 10 
14 


0.40 0.165 0.019 0.009 0.002 1.13 
0.28 0.280 0.020 0.006 0.003 1.05 


-10 10 

14 


0.295 0.096 018 0.007 0.002 0.48 
0.220 0.130 0.019 0.007 0.003 0.44 


oo 10 
14 


0.315 0.103 0.020 0.010 0.003 0.50 
0.181 0.102 0.024 0.006 0.004 0.44 



TABLE II: Values of kl in units of 



akp 


N=ll,13 


N=15,17,19 





1 


2 


-1 


1 


2 


-3 


1 


2 


-10 


1 


1 


oo 


1 


1 


10 


1 


1 


3 


1 


1 


2 


or 1 


or 1 


1 








0.75 








0.5 
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TABLE III: Summary of the results in the strongly interacting regime. 



1 

akp 


A 

EpG 


A EpG 


2 EpG 


R-rms 

ro 




U 


n no f A\ 


U. 44(^1 j 


U 


oo 


OO 


n 1 
U.i 


i.Uo(^0 ) 


U.o4(^i j 




o.oy 


44. o 


0.3 


1.37(5) 


0.02(1) 


-0.20(1) 


1.21 


14.5 


0.5 


1.80(5) 


-0.33(1) 


-0.49(1) 


0.74 


8.9 


1.0 


3.2(1) 


-2.23(1) 


-2.31(1) 


0.38 


4.6 


1.3 


5.7(3) 


-4.58(2) 


-4.63(1) 


0.28 


3.4 


2.0 


14.0(5) 


-12.84(3) 


-12.86(1) 


0.19 


2.3 
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